#RScript
#Local Competitive Authoritarianism and Post-Conflict Violence.An Analysis of the Assassination of Social Leaders in Colombia 
#Juan Albarracin, Juan Pablo Milanese, Inge Valencia, and Jonas Wolff
#July 28, 2022


##Set working directory to store results.
##Load data: 
library(haven)
lideres <- read_dta("lideres.dta")
 
#Load packages, if not installed, use install.packages() 
library(MASS)
library(reshape2) 
library(ggplot2)
library(nnet)
library(psych)
library(outreg)

##Preparing data

#Log Transformations
##Coca Hectares
lideres$logcoca<-log(lideres$hectaresas_coca16 + 1)
##Gold production (in ounzes)
lideres$gold<-log(lideres$onzasoro17 + 1)


########### TABLE 1: Descriptive Statistics ###############

descriptives<-describe(lideres[ , c("asesinatos_total", "asesinatos_new", "logleft2015", "Participacion", "nep_mun", "integral", "solicitud_RTDAF", "FARC_presence", "multiple_presence", "logcoca", "gold", "tasahom1000", "log_poblacion")],fast=TRUE)
write.csv(descriptives, "table1.csv")

#Transform variables to factor variables
lideres$FARC_presence<-as.factor(lideres$FARC_presence)
lideres$Guerrilla_presence<-as.factor(lideres$Guerrilla_presence)
lideres$paras_presence<-as.factor(lideres$paras_presence)
lideres$criminal_presence<-as.factor(lideres$criminal_presence)

##########################TABLE 2##################################################

##################Logit Models ####################

##Table 2, Model 1.

###Create object with regression results
resultados_logit<-glm(asesinatos_new ~ logleft2015 + Participacion + nep_mun + integral + solicitud_RTDAF + FARC_presence + multiple_presence +  logcoca + gold + tasahom1000 + log_poblacion,
                       family=binomial, data=lideres)
###View Results
summary(resultados_logit)

##Table 2, Model 2. 
###Create object with regression results
resultados_logit_quadratic<-glm(asesinatos_new ~  logleft2015 + I(logleft2015^2) + Participacion + nep_mun + I(nep_mun^2)  +  integral + solicitud_RTDAF + FARC_presence + multiple_presence +  logcoca + gold + tasahom1000 + log_poblacion,
                                family=binomial, data=lideres)
###View Results
summary(resultados_logit_quadratic)

############### Count Models ######################

#Model 3, Table 2
###Create object with regression results
count <- glm.nb(asesinatos_total ~ logleft2015 + Participacion + nep_mun +  integral + solicitud_RTDAF + FARC_presence + multiple_presence +  logcoca + gold + tasahom1000 + log_poblacion, data=lideres)
###View results
summary(count)

#Model 4, Table 2
###Create object with regression results
count_quadratic <- glm.nb(asesinatos_total ~ logleft2015 + I(logleft2015^2) + Participacion + nep_mun + I(nep_mun^2)  +  integral + solicitud_RTDAF + FARC_presence + multiple_presence +  logcoca + gold + tasahom1000 + log_poblacion, data=lideres)
###View results
summary(count_quadratic)

##################Create object with table regression results################################3

table2<-outreg(list("Model 1 Logit" =resultados_logit, "Model 2 Logit"=resultados_logit_quadratic, "Model 3 Negtive Binomial"=count, "Model 4 Negative Binomial"=count_quadratic), digits = 4L, alpha = c(0.05, 0.01, 0.001), 
               bracket = c("se"), starred = c("coef"), robust = FALSE, small = TRUE,
               constlast = FALSE, norepeat = TRUE)
write.csv(table2, "table2.csv")


############ GRAPHS based on Model 2, Table 2 (logistic model) ####################

## Figure 1: Share of Leftist Vote and the Probability of an Assassination
###Create object with data to graph
newdata7 <- data.frame(nep_mun = mean(lideres$nep_mun, na.rm = TRUE), Participacion = mean(lideres$Participacion, na.rm = TRUE), logleft2015 = rep(seq(0,5,0.5)), integral = mean(lideres$integral, na.rm = TRUE), solicitud_RTDAF=mean(lideres$solicitud_RTDAF, na.rm = TRUE), FARC_presence="1", multiple_presence=1, logcoca=mean(lideres$logcoca, na.rm = TRUE), gold=mean(lideres$gold, na.rm = TRUE), tasahom1000=mean(lideres$tasahom1000, na.rm = TRUE), log_poblacion=mean(lideres$log_poblacion, na.rm = TRUE))
newdata7 <- cbind(newdata7, predict(resultados_logit_quadratic, newdata = newdata7, type = "link",
                                    se = TRUE))
newdata7 <- within(newdata7, {
  PredictedProb <- plogis(fit)
  LL <- plogis(fit - (1.96 * se.fit))
  UL <- plogis(fit + (1.96 * se.fit))
})
###Plot using ggplot
ggplot(newdata7, aes(x = logleft2015, y = PredictedProb)) + geom_ribbon(aes(ymin = LL, ymax = UL), alpha = 0.2) + geom_line(size = 1) + labs(x = "% vote for the left (log)", y = "Probability of at least one assassination") + theme_bw()
ggsave("figure1.jpeg", dpi = 600 )

## Figure 2: Effective Number of Parties and Probability of an Assassination
###Create object with data to graph
newdata5 <- data.frame(nep_mun = rep(seq(1,11,0.5),2), Participacion = mean(lideres$Participacion, na.rm = TRUE), logleft2015 = mean(lideres$logleft2015,na.rm = TRUE), integral = mean(lideres$integral, na.rm = TRUE), solicitud_RTDAF=mean(lideres$solicitud_RTDAF, na.rm = TRUE), FARC_presence="1", multiple_presence=1, logcoca=mean(lideres$logcoca, na.rm = TRUE), gold=mean(lideres$gold, na.rm = TRUE), tasahom1000=mean(lideres$tasahom1000, na.rm = TRUE), log_poblacion=mean(lideres$log_poblacion, na.rm = TRUE))
newdata5 <- cbind(newdata5, predict(resultados_logit_quadratic, newdata = newdata5, type = "link",
                                    se = TRUE))
newdata5 <- within(newdata5, {
  PredictedProb <- plogis(fit)
  LL <- plogis(fit - (1.96 * se.fit))
  UL <- plogis(fit + (1.96 * se.fit))
})
###Plot using ggplot
ggplot(newdata5, aes(x = nep_mun, y = PredictedProb)) + geom_ribbon(aes(ymin = LL,
                                                                        ymax = UL), alpha = 0.2) + geom_line(size = 1) + labs(x = "Effective number of parties", y = "Probability of at least one assassination") + theme_bw()

ggsave("figure2.jpeg", dpi = 600 )

##Figure 3: Turnout and the Probability of an Assassination
###Create object with data to graph
newdata8 <- data.frame(nep_mun = mean(lideres$nep_mun, na.rm = TRUE), Participacion = rep(seq(20,100,5),2), logleft2015 = mean(lideres$logleft2015,na.rm = TRUE), integral = mean(lideres$integral, na.rm = TRUE), solicitud_RTDAF=mean(lideres$solicitud_RTDAF, na.rm = TRUE), FARC_presence="1", multiple_presence=1, logcoca=mean(lideres$logcoca, na.rm = TRUE), gold=mean(lideres$gold, na.rm = TRUE), tasahom1000=mean(lideres$tasahom1000, na.rm = TRUE), log_poblacion=mean(lideres$log_poblacion, na.rm = TRUE))
newdata8 <- cbind(newdata8, predict(resultados_logit_quadratic, newdata = newdata8, type = "link",
                                    se = TRUE))
newdata8 <- within(newdata8, {
  PredictedProb <- plogis(fit)
  LL <- plogis(fit - (1.96 * se.fit))
  UL <- plogis(fit + (1.96 * se.fit))
})
###Plot using ggplot
ggplot(newdata8, aes(x = Participacion, y = PredictedProb)) + geom_ribbon(aes(ymin = LL,
                                                                              ymax = UL), alpha = 0.2) + geom_line(size = 1) + labs(x = "Turnout (%)", y = "Probability of at least one assassination") + theme_bw()

ggsave("figure3.jpeg", dpi = 600 )


############ GRAPHS based on Model 4, Table 2 (count model) ####################

## Figure 4: Share of Leftist Vote and Predicted Assassination Count
###Create object with data to graph
newdata3 <- data.frame(nep_mun = mean(lideres$nep_mun, na.rm = TRUE), Participacion = mean(lideres$Participacion, na.rm = TRUE), logleft2015 = seq(0,5,0.5), integral = mean(lideres$integral, na.rm = TRUE), solicitud_RTDAF=mean(lideres$solicitud_RTDAF, na.rm = TRUE), FARC_presence="1", multiple_presence=1, logcoca=mean(lideres$logcoca, na.rm = TRUE), gold=mean(lideres$gold, na.rm = TRUE), tasahom1000=mean(lideres$tasahom1000, na.rm = TRUE), log_poblacion=mean(lideres$log_poblacion, na.rm = TRUE))
newdata3 <- cbind(newdata3, predict(count_quadratic, newdata3, type = "link", se.fit=TRUE))
newdata3 <- within(newdata3, {
  Assassination <- exp(fit)
  LL <- exp(fit - (1.96 * se.fit))
  UL <- exp(fit + (1.96 * se.fit))
})
###Create object with data to graph
ggplot(newdata3, aes(logleft2015, Assassination)) +
  geom_ribbon(aes(ymin = LL, ymax = UL),  alpha = .25) +
  geom_line(size = 1) +
  labs(x = "% vote left (log)", y = "Predicted number of assassinations") + theme_bw()
ggsave("figure4.jpeg", dpi = 600 )

## Figure 5: Effective Number of Parties and Predicted Assassination Count
###Create object with data to graph
newdata4 <- data.frame(nep_mun = seq(1,10,0.5), Participacion = mean(lideres$Participacion, na.rm = TRUE), logleft2015 = mean(lideres$logleft2015, na.rm = TRUE), integral = mean(lideres$integral, na.rm = TRUE), solicitud_RTDAF=mean(lideres$solicitud_RTDAF, na.rm = TRUE), FARC_presence="1", multiple_presence=1, logcoca=mean(lideres$logcoca, na.rm = TRUE), gold=mean(lideres$gold, na.rm = TRUE), tasahom1000=mean(lideres$tasahom1000, na.rm = TRUE), log_poblacion=mean(lideres$log_poblacion, na.rm = TRUE))

newdata4 <- cbind(newdata4, predict(count_quadratic, newdata4, type = "link", se.fit=TRUE))
newdata4 <- within(newdata4, {
  Assassination <- exp(fit)
  LL <- exp(fit - (1.96 * se.fit))
  UL <- exp(fit + (1.96 * se.fit))
})
###Create object with data to graph
ggplot(newdata4, aes(nep_mun, Assassination)) +
  geom_ribbon(aes(ymin = LL, ymax = UL), alpha = .25) +
  geom_line(size = 1) +
  labs(x = "Effective number of parties", y = "Predicted number of assassinations") + theme_bw()
ggsave("figure5.jpeg", dpi = 600 )

## Figure 6: Turnout and Predicted Assassination Count
###Create object with data to graph
newdata2 <- data.frame(nep_mun = mean(lideres$nep_mun, na.rm = TRUE), Participacion = seq(20,100,5), logleft2015 = mean(lideres$logleft2015, na.rm = TRUE), integral = mean(lideres$integral, na.rm = TRUE), solicitud_RTDAF=mean(lideres$solicitud_RTDAF, na.rm = TRUE), FARC_presence="1", multiple_presence=1, logcoca=mean(lideres$logcoca, na.rm = TRUE), gold=mean(lideres$gold, na.rm = TRUE), tasahom1000=mean(lideres$tasahom1000, na.rm = TRUE), log_poblacion=mean(lideres$log_poblacion, na.rm = TRUE))
newdata2 <- cbind(newdata2, predict(count_quadratic, newdata2, type = "link", se.fit=TRUE))
newdata2 <- within(newdata2, {
  Assassination <- exp(fit)
  LL <- exp(fit - (1.96 * se.fit))
  UL <- exp(fit + (1.96 * se.fit))
})
###Create object with data to graph
ggplot(newdata2, aes(Participacion, Assassination)) +
  geom_ribbon(aes(ymin = LL, ymax = UL), alpha = .25) +
  geom_line(size = 1) +
  labs(x = "% turnout", y = "Predicted number of assassinations") + theme_bw()

ggsave("figure6.jpeg", dpi = 600 )

########################  APPENDIX  ####################################
#############Additional Analyses for the Appendix ######################

#########################  TABLE A1 ######################################

#Model A1: Different Measure of Competitiveness
resultados_logit_diff<-glm(asesinatos_new ~ logleft2015 + I(logleft2015^2) +  Participacion + DifCon + integral + solicitud_RTDAF + FARC_presence + multiple_presence +  logcoca + gold + tasahom1000 + log_poblacion,
                            family=binomial, data=lideres)
summary(resultados_logit_diff)

#Model A2: With Dummy Variables for Other Armed Actors
resultados_logit_otros<-glm(asesinatos_new ~  logleft2015 + I(logleft2015^2) + Participacion + nep_mun + I(nep_mun^2)  +  integral + solicitud_RTDAF + FARC_presence + Guerrilla_presence + paras_presence + criminal_presence + logcoca + gold + tasahom1000 + log_poblacion,
                            family=binomial, data=lideres)

##Model A3: With Local Vanhanen as an Alternative Measure of Local Competitive Authoritarianism
resultados_logit_vanhanen<-glm(asesinatos_new ~  logleft2015 + Vanhanen + integral + solicitud_RTDAF + FARC_presence + multiple_presence +  logcoca + gold + tasahom1000 + log_poblacion,
                               family=binomial, data=lideres)
summary(resultados_logit_vanhanen)

##Model A4: With Local Vanhanen (Quadratic) as an Alternative Measure of Local Competitive Authoritarianism
resultados_logit_vanhanen2<-glm(asesinatos_new ~  logleft2015 + I(logleft2015^2) + Vanhanen + I(Vanhanen^2) +  integral + solicitud_RTDAF + FARC_presence + multiple_presence +  logcoca + gold + tasahom1000 + log_poblacion,
                                family=binomial, data=lideres)
summary(resultados_logit_vanhanen2)

##Create object with table results for all models in A1
tableA1<-outreg(list(resultados_logit_diff, resultados_logit_otros, resultados_logit_vanhanen, resultados_logit_vanhanen2), digits = 3L, alpha = c(0.05, 0.01, 0.001), 
               bracket = c("se"), starred = c("coef"), robust = FALSE, small = TRUE,
               constlast = FALSE, norepeat = TRUE)
write.csv(tableA1, "tableA1.csv")


#########################  TABLE A2 ######################################

#Models with different measurement of state capacity (judicial effectiveness)

##Model A5, Table A2: Logit
resultados_logit_basic_extra<-glm(asesinatos_new ~  logleft2015 + Participacion + nep_mun +  judicial_mean + solicitud_RTDAF + FARC_presence + multiple_presence +  logcoca + gold + tasahom1000 + log_poblacion,
                                  family=binomial, data=lideres)
summary(resultados_logit_basic_extra)

##Model A6, Table A2: Logit with quadratic terms
resultados_logit_quadratic_extra<-glm(asesinatos_new ~  logleft2015 + I(logleft2015^2) + Participacion + nep_mun + I(nep_mun^2)  +  judicial_mean + solicitud_RTDAF + FARC_presence + multiple_presence +  logcoca + gold + tasahom1000 + log_poblacion,
                                      family=binomial, data=lideres)
summary(resultados_logit_quadratic_extra)

##Model A7, Table A2: Count model
count_extra <- glm.nb(asesinatos_total ~ logleft2015 + Participacion + nep_mun +  judicial_mean + solicitud_RTDAF + FARC_presence + multiple_presence +  logcoca + gold + tasahom1000 + log_poblacion, data=lideres)
summary(count_extra)

##Model A8, Table A2: Count model with quadratic terms
count_quadratic_extra <- glm.nb(asesinatos_total ~ logleft2015 + I(logleft2015^2) + Participacion + nep_mun + I(nep_mun^2)  +  judicial_mean + solicitud_RTDAF + FARC_presence + multiple_presence +  logcoca + gold + tasahom1000 + log_poblacion, data=lideres)

##Create object with table results for all models in A2
tableA2<-outreg(list(resultados_logit_basic_extra, resultados_logit_quadratic_extra, count_extra, count_quadratic_extra), digits = 3L, alpha = c(0.05, 0.01, 0.001), 
                bracket = c("se"), starred = c("coef"), robust = FALSE, small = TRUE,
                constlast = FALSE, norepeat = TRUE)
write.csv(tableA2, "tableA2.csv")

#########################  TABLE A3 ######################################
################Results for assasinations ocurring in years 2017-2019 #####

##Model A9, Table A3
resultados_logit_as2<-glm(asesinatos_dummy_1719 ~ logleft2015 + Participacion + nep_mun +  integral + solicitud_RTDAF + FARC_presence + multiple_presence +  logcoca + gold + tasahom1000 + log_poblacion,
                          family=binomial, data=lideres)
summary(resultados_logit_as2)

##Model A10, Table A3
resultados_logit_quadratic2<-glm(asesinatos_dummy_1719 ~ logleft2015 + I(logleft2015^2) + Participacion + nep_mun + I(nep_mun^2)  +  integral + solicitud_RTDAF + FARC_presence + multiple_presence +  logcoca + gold + tasahom1000 + log_poblacion,
                                 family=binomial, data=lideres)
summary(resultados_logit_quadratic2)

##Model A11, Table A3
count2 <- glm.nb(asesinatos_total_1719 ~ logleft2015 +  Participacion + nep_mun + integral + solicitud_RTDAF + FARC_presence + multiple_presence +  logcoca + gold + tasahom1000 + log_poblacion, data=lideres)
summary(count2)

##Model A12, Table A3
count_quadratic2 <- glm.nb(asesinatos_total_1719 ~ logleft2015 + I(logleft2015^2) + Participacion + nep_mun + I(nep_mun^2)  +  integral + solicitud_RTDAF + FARC_presence + multiple_presence +  logcoca + gold + tasahom1000 + log_poblacion, data=lideres)
summary(count_quadratic2)

##Create object with table results for all models in A3
tableA3<-outreg(list(resultados_logit_as2, resultados_logit_quadratic2, count2, count_quadratic2), digits = 4L, alpha = c(0.05, 0.01, 0.001), 
                bracket = c("se"), starred = c("coef"), robust = FALSE, small = TRUE,
                constlast = FALSE, norepeat = TRUE)
write.csv(tableA3, "tableA3.csv")
